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Abstract 

The finite basis optimized effective potential (OEP) method within density functional theory is 
examined as an ill-posed problem. It is shown that the generation of nonphysical potentials is a 
controllable manifestation of the use of unbalanced, and thus unsuitable, basis sets. A modified 
functional incorporating a regularizing smoothness measure of the OEP is introduced. This pro- 
vides a condition on balanced basis sets for the potential, as well as a method to determine the 
most appropriate OEP potential and energy from calculations performed with any finite basis set. 

PACS numbers: 31.10,+z, 31.15Ew, 31.15.Pf 
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Kohn-Sham density functional theory (KS DFT) [l] enjoys wide application owing to 
the computational accessibility afforded by formulating the many body problem in terms 
of the noninteracting KS reference system. However, the exact exchange-correlation energy 
functional E xc and corresponding local potential, v xc , are unknown. The future success of 
DFT is dependent on the availability of suitable approximations for E xc . Significant interest 
is being shown for the development of implicit density functionals, depending explicitly 
on the KS orbitals [2]. For such functionals, v xc cannot be directly obtained as a simple 
functional derivative and requires an OEP method for its determination 0,3, 51. 

The OEP concept first appeared within the Hartree-Fock (HF) formalism fjj, and later 
was employed in DFT [7] using the HF exact exchange energy functional (EXX). These 
conventional approaches identify the local KS potential for exact exchange as that obtained 
from solving 5E[{<fif s }]/5v s (r) = 0, which in turn leads to a linear integral equation to be 
solved for the OEP. The rigorous justification of the OEP method as a variational minimiza- 
tion has however only recently been provided within the potential functional formulation 
8j. This formulation lends itself to a direct approach to solving for the OEP 5] where the 
minimization of the energy functional over local potentials is considered. 

The discrete representation of the OEP problem in finite basis sets for both the Kohn- 
Sham orbitals and the potential can be ill-posed 0, [hJ: While the total OEP energy is 
stable with respect to the changes in the potential, there can be many different potentials, 
including nonphysical ones, having numerically degenerate total energies. This has lead to 
some degree of confusion of what in fact constitutes a valid, finite basis, OEP implementation 
[h]], when the ill-posedness is not appropriately accounted for. 

In this Letter, we will show that the ill-posed nature of the discrete OEP originates from 
the use of basis sets that are unsuitable, as a result of being unbalanced. This ill-posedness 
does not imply that the OEP method is unphysical by construction, rather, just that it can 
lead to nonphysical potentials. The regularization method developed in this work insures 
that the physical context of finite basis OEP calculations are maintained in all cases, with 
generation of physically meaningful potentials. 

The following analysis reveals the physical origin of the ill-posed nature in the discrete 
OEP problem. In determining the OEP, one looks for the total-energy minimizing po- 
tential v s (r) whose lowest N eigenstates are the set of wavefunctions with the finite basis 
set expansion, \<pi) = c «mIXm)' f° r i — 1j •••> satisfying the one-electron equations 
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Y.^(Xv\ - 7>V 2 + v s (r)\xn)ciii = ^iJ2^(Xu\Xn)ci fi - This equation is the only link from v s (r) 
to the OEP energy. If a variation of v s (r) is outside of the range of the orbital basis set 
{X/i} ~ meaning such a variation in the potential leads to numerically the same total energy 
and orbitals \<pi) - then this variation can lead to a nonphysical potential. The nonphysical 
potentials reported in Ref. [h| appear to come from such a construction. For nonphysi- 
cal potentials with wild oscillations, the orbital basis set is certainly inadequate for solving 
the one-electron equations and the discrete OEP equation is a poor approximation to the 
original OEP problem. Given an orbital basis set {Xfi}i there are many variations of v s (r) 
that can lead to such behavior, one example being adding to v s (r) Gaussian functions with 
exponents much higher than the exponents in the orbital basis set {Xy}- We consider such 
a potential basis set unbalanced with the orbital basis set. 

While a careful choice of basis sets for the orbitals and the potential can lead to accurate 
potentials, as shown in many previous calculations 4], [l], El, we aim here to develop a robust 
method which can deliver accurate OEP potentials and energies, with both balanced and 
unbalanced basis sets. 

To examine the reliability of our approach we use the OEP procedure with the LDA 
(SVWN5) functional so that a direct comparison between the calculated OEP potentials 
and the LDA potential can be made. The success of our method for LDA OEP transfers to 
EXX as will be shown. 

Our OEP implementation is the direct optimization approach of Yang and Wu [s| where 
the trial potential is expanded in a finite basis set, {gt}, as v° (r) = v ext (r)+vo(r)+J2t ^t9t( T )- 
Here i> ex t is the external potential of the system under consideration and vq is a fixed reference 
potential, taken as the Fermi- Amaldi potential (or the Coulomb potential for LDA) for the 
sum of the atomic densities so to enforce the correct asymptotic behavior upon v s . Transfer- 
ring the functional dependence from the KS potential on to the expansion coefficients {b°} 
in this way gives rise to an efficient implementation of the OEP based on the unconstrained 
minimization of E({b1}) with readily available analytic derivatives The iterative mini- 



mization of E({b%}) is efficiently achieved using an approximate Newton method rjj which 
takes an approximate form for the Hessian as H£ t = J2 ia ^ og|g " l ^^ ag|gt| ^ g> . 

In the following we will consider two forms of regularization for extracting physical po- 



tentials from ill-posed finite basis OEP calculations. The first, used by us before [5|, 111] , 
however will be shown to be inadequate to remove all irregularities in the resulting potential. 
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We then introduce a new form of regularization for the OEP that produces potentials of the 
highest quality, along with the corresponding energy. 

We will focus on the specific example of argon with the orbital basis cc-pVDZ and consider 
the effects of different potential basis on the potentials and energies obtained from the 
OEP procedure. The specific basis sets we use for the potential are the orbital basis itself, 
and three s-type even tempered (ET) basis sets denoted as ArX (X=64, 1024, 8192) with 
exponents 2 n , —4 < n < log 2 (X). 
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FIG. 1: Variation in the smoothness of the LDA OEP potential as a function of the SVD cutoff 
for argon (orbital basis cc-pVDZ) with four different potential basis sets. Insert: Spectrum of the 
approximate Hessian for each basis set. 

The nature of the Hessians spectrum is intimately related to the behavior of the potentials 
obtained. Most applications of the direct OEP procedure have made use of coinciding 
orbital-quality basis sets for both orbitals and potential. This, generally, generates physical 
potentials and energies. From the Hessian spectrum in Fig. 1 we see that using the cc-pVDZ 
orbital and potential basis leads to a full rank Hessian resulting in the convergence of the 
OEP to a physical potential. However, if we use other constructions for the potential basis 
such as ET basis sets, we can observe much different spectral structures for the Hessian. 
The rank deficiency for the three ET basis sets is obvious from the significant gap in the 
spectra of Fig. 1. To solve the update equation in an stable manner, both truncated 
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FIG. 2: TSVD-regularized LDA potentials of argon for various values of the SVD cutoff, a, using 
the potential basis Ar8192. The solid line depicts the exact LDA potential for comparison. 



SVD (TS VD) and Tikhonov regularizations have been proposed and implemented llj . The 
TSVD regularization, to be used in this discussion, neglects contributions from eigenvectors 
corresponding to those singular values below some chosen cutoff. Fig. 1 shows how the 
smoothness of the resulting OEP, as measured by the norm ||V^ t 6tfl't( r )||; varies as a 
function of the TSVD cutoff used when solving the Hessian update equation. As the cutoff 
approaches the lower cluster of singular values, the structure becomes step-like, with jumps 
in the norm as the cutoff passes through eigenvalues of the Hessian. The potential becomes 
very nonphysical (Fig. 2), reminiscent of those displayed in Ref. [h]]. It is clear that this 
TSVD regularization, on its own, is not sufficient to ensure a physically meaningful solution 
for the OEP in general. It does however confine irregularities to regions near the nuclei, 
which could be sufficient in many applications. 



Next, we take the classical approach to il 
desirable measure to regularize the solution 



posed problems and further incorporate some 



121 ] . It is clear that any nonphysical oscillatory 



behavior in the potential will be confined to the basis set expansion t>b = Y^t^tdt- We 
thus introduce a A-regularization by constraining our solutions to yield smooth potentials 
as measured by the smoothing norm ||Vt>b||, thereby restricting the nonphysical variations 
in the potential because of unbalanced basis set. This norm is certainly not unique, however 
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is simple to implement and will be seen to produce very satisfying results. We define a 
regularized energy functional as 

^ A (b) = i^ YW (b) + A||V^ b (r)|| 2 , (1) 

where ||Vt>b(r)|| 2 = 2b T Tb, T is the kinetic energy integral matrix in the potential basis 
and E YW (h) is the OEP energy calculated according to the Yang-Wu method. The energy 
derivatives with respect to the coefficients are modified accordingly as Vb^A = Vb-E YW + 
4ATb and V b ^A = V b -E YW + 4AT. The significance of the smoothing contribution relative 
to the energy minimization is then determined by the magnitude of A. Our method of 
determining the most appropriate value of A (that is, the value of A which corresponds to a 
physical potential) is the classical L-curve analysis [l2j. This involves constructing the plot 
of log(||Vfb||) against log(E YW — E rei ) for A over some appropriate interval. The reference 
energy, £ ref , is taken as E UF for EXX calculations and E LBA for LDA OEP. 

This curve presents the trade-off between the desire to extract a smooth potential, with 
that of an unconstrained minimization of the energy functional which can allow the intro- 
duction of nonphysical potentials. For a balanced potential basis set - defined here as one 
that converges to a smooth physical potential - gradually increasing A from zero has no 
significant effect on the norm ||Vi>b||- In such a situation no A-regularization is needed. On 
the other hand, a general basis can see a significant decrease in the norm as A increases, 
until there comes a region where variations in the the norm subside corresponding to the 
formation of a stable and physical potential. 

The most appropriate choice of A for any particular physical application is not always a 
trivial decision, dependent on the particular system under consideration and the tolerance 
that are acceptable [l2j|. In the present context, our calculations have lead to a consistent 
definition: The A corresponding to the optimal potential, and thus also the OEP energy for 
the given finite basis sets and chosen smoothing norm, is that point on the L-curve with 
minimum slope and, if it is not unique, that has the minimum value of ||Vt>b||- Support for 
such a choice is presented in 

Figs. 3 and 5 present such L-curve analysis for our argon example with LDA and EXX 
functionals. We use the notation AO/PB(A) to describe calculations performed with AO and 
PB atomic orbital and potential basis sets respectively, and with a regularization parameter 
of A. For consistency we chose to perform the calculations without any TSVD regularization. 
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FIG. 3: Argon (cc-pVDZ) LDA OEP L-curve for various potential basis sets. 

The curve for cc-pVDZ/cc-pVDZ (A) shows that no regularization needs to be applied to 
solve for an appropriate OEP as explained previously, however the cost of using such a basis 
is an energy higher than that can be obtained from more flexible ET basis sets. The most 
appropriate LDA OEP potential for this basis set (cc-pVDZ/cc-pVDZ (0)) as seen in Fig. 
4 underestimates the LDA potential near the nucleus. 

All three of the ET basis sets produce poor potentials with oscillations extending far from 
the nuclei when no regularization is used. The LDA L-curves for cc-pVDZ/Ar8192(A) and 
cc-pVDZ/Arl024(A) are of an ideal form - as A tends towards zero smooth potentials are 
generated with numerical energy agreement with the LDA value, at which point the curve 
rises vertically corresponding to the introduction of nonphysical oscillations in the potential. 
In this situation the choice of a meaningful potential is clear, the corner of the L-curve is 
well defined and A* = 10~ 14 , the optimal choice of A, is that point infinitesimally close to the 
corner (with minimum slope). The corresponding potential (Fig. 4) agrees with the LDA 
potential in all but a small region around the nuclei. 

The EXX results for the same basis sets in Figs. 5 and 6 show that the critical value 
of A occurs at ~ AE = 10 -4 a.u. Significantly, this example illustrates that nonphysical 
potentials can be obtained when the ill-posedness is not addressed, and that the physical 
OEP potentials can be obtained with our method. As another example, the EXX OEP 
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FIG. 4: Optimal LDA OEP potentials for argon (cc-pVDZ) obtained from the L-curve analysis for 
each basis set. 
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FIG. 5: Argon (cc-pVDZ) EXX L-curve for various potential basis sets. 



L-curve analysis for N 2 is shown in Fig. 7. We see that all basis sets regularize remarkably 
consistently, giving rise to a consistent OEP energy that will be explored in some detail 
within a future paper. Further detailed examples for atoms and molecules are provided 
online via EPAPS [13 1. 
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FIG. 6: Optimal EXX OEP potentials for argon (cc-pVDZ) obtained from the L-curve analysis for 
each basis set. 

In conclusion, we have constructed a regularization procedure by introducing a measure 
of the smoothness of the OEP into the energy functional. This, together with an L-curve 
analysis, allows us to determine the physically meaningful OEP potential and energy from 
calculations performed with any finite basis sets. We have shown that one quantitative 
measure of such balance is provided by the approximate Hessian, however this does not 
provide an a priori judgement needed for routine application of finite basis OEP calculations. 
For this, the construction of balanced potential basis sets are required. Support from the 
National Science Foundation is gratefully acknowledged. 
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Supplementary Material 



I. OUR CHOICE OF THE OPTIMAL REGULARIZATION PARAMETER A* 



Firstly, we wish to consider our choice of A*: The value of the regularization parame- 
ter corresponding to the optimal choice of potential. In many statistical applications, it 
is recommended to choose that corner of the L-curve with maximum curvature, for this 
point is argued to represent the optimal tradeoff between stability of solutions, and absolute 
minimization of the quantity of interest 12]. This however is dependent on the physical 
system under consideration and the tolerances acceptable to the user. We find the maxi- 
mum curvature point on the L-curve not to be optimal for our concern - irregularities have 
been introduced and the regularized potential has began to diverge away from the physical 
potential. The over-smoothing at the corner of the L curve has also been observed in some 
applications [ijj]. 

The following example shows, for an LDA calculation of nitrogen (potential basis set is 
introduced in section III), the changes in the potential as we traverse the corner of the L-curve 
in Fig. (IE]). A = 10~ 4 best corresponds to that point on the L-curve with minimum slope, 
(or minimum absolute gradient), while the maximum curvature corresponds to A = 10~ 6 . 
The irregularities appearing with the different A are most noticeable at the origin and for 
large z. 

In the following we will go though several more examples showing the versatility of our 
procedure. 



II. WATER 



We shall consider the orbital basis of cc-pVDZ, and the three potential basis sets cc-pVDZ, 
H2O01, and H2O02. H2O01 and H2O02 are even tempered (ET) basis sets constructed as 
follows: H2O01 - 9s9p (O) and 9s (H) uncontracted cartesian Gaussians with exponents 
2 n , —3 < n < 5. H2O02 - 18sl5p (O) and 18s (H) uncontracted cartesian Gaussians with 
exponents 2™, —4 < n < 13 (s) and 2 n , —4 < n < 13 (p). Our construction of the ET basis 
sets are such to encourage poor behavior in the resulting, unregularized, potentials so to 
show the versatility of our procedure. 
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FIG. 8: LDA OEP potentials for N2 (cc-pVDZ) in the vicinity of the L-curve corner for the potential 
basis N2256. 

We use the geometry 0.917A for the O-H bond length and 104.4 degrees for the HOH 
angle, and plot potentials along the HOH bisector. 

A. LDA 

For water, the Hessian spectra are shown in Fig. (Q. As is typical for molecules, it 
becomes much harder to interpret the Hessian spectra in terms of a definable cluster of 
singular values that we would consider to be associated with null-space eigenvectors. The 
problem becomes very much ill-conditioned for the basis H2O02: the eigenvectors decay 
gradually to (numerically) zero, as opposed to the situation presented in the paper for Ar 
where the problem is very much rank deficient. 

The variation of 1 1 Vi>b| | as a function off the TSVD cutoff, a*, is presented in Fig. (jTOl) . in 
an attempt for extract a definable range of eigenvalues for which we can consider the Hessian 
to be stable if used in solving the update equation H(b n )p = —VE(b n ). The step structure 
as was seen in the main text is still present, however for the H2O02 basis it becomes very 
hard to argue that any particular value of the TSVD cutoff would be any more favorable for 
a TSVD-only regularized solution for the OEP. The other two basis sets provide very clear 
interpretations in terms of the eigenspectrum of Fig. ([9]). H2O01 has a smallest eigenvalue 
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FIG. 9: Spectrum of the approximate Hessian for each basis set. 



of only 3.5 x 10~ 12 , so that any choice for a TSVD cutoff smaller than this will yield identical 
results. The resulting potential is, however, very poor with nonphysical oscillations about 
the oxygen atom as seen in Fig. ffTTj) . Choosing a value for the cutoff around 10~ 5 does lead 
to much better potentials, but still there are some irregularities around the nucleus. The 
situation for cc-pVDZ/cc-pVDZ is typical, the Hessian eigenvalues decay to only ~ 3 x 10~ 5 
a.u. 

It is clearly necessary to apply the A-regularization for both ET basis sets in order to 
produce appropriate potentials. We chose not to apply any TSVD when solving the up- 
date equation and obtain the L-curves seen in Fig. (fT2|) . in turn giving the potentials, 
corresponding to the minimum slope point, in Fig. ( TT3j) 



The Hessian structure is very similar to that of the LDA case. This is a very general 
observation for all the systems we have examined, so the discussion for the LDA case carries 
over to EXX. As such we will specifically only discuss the L-curves and the potentials we 
can extract in the EXX case although a complete set of figures are presented. 

The L-curves of Fig. f|T7|) agree exceptionally well for all but the smallest amounts of 



B. EXX 
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FIG. 11: TSVD-regularized LDA potentials for water (cc-pVDZ) for various values of the cutoff. 
The solid line depicts the exact LDA potential. 

regularization. This is very encouraging given the differences in the potential basis sets. An 
inflexion point is well defined and gives an approximation to the OEP energy of E OEP — 
E RF ~ 1 x 10~ 25 a.u.. If we want to get more accurate energies a finer A-sampling grid would 
be necessary. The given basis sets can reproduce the HF energy to an agreement of ~ 10~ 5 
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FIG. 13: Optimal LDA OEP potentials for water (cc-pVDZ) obtained from the L-curve analysis 
for each basis set. 

a.u. with nonphysical potentials. The potentials we extract from the L-curve are shown in 
Fig. (TT5|) . The two ET basis sets produce essentially identical potentials, differing only at 
the oxygen cusp. The potential cc-pVDZ/cc-pVDZ(0) is also of an acceptable quality, but 
as expected with the use of basis sets designed for orbitals, it does struggle to reproduce the 
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FIG. 14: Spectrum of the approximate Hessian for each basis set. 
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finer details that a flexible ET basis can. 



16 



■cc-pVDZ/H2OO1(0*=O) 
■cc-pVDZ/H2O01(o*=10"°) 
cc-pVDZ/H2O01(o*=10" 4 ) 




FIG. 16: TSVD-regularized EXX potentials for water (cc-pVDZ/H2O01) for various values of the 
TSVD cutoff 



3.5 
3 

2.5 
2 

5> 1.5 



g 1 1 



0.5 


-0.5 



A □ 
An 



□ cc-pVDZ/H2O01(Ji) 
A cc-pVDZ/H2O02(^) 



-4.5 -4 -3.5 -3 -2.5 
log(E YW -^ F ) 



-2 -1.5 -1 



FIG. 17: LDA L-curves for water in the various potential basis sets. 



III. NITROGEN 



The potential basis sets we use for N2 are as denoted as N2256, N2512, and N28192, 
consisting of 13s3p, 14s5p and 18s8p uncontracted cartesian Gaussians with exponents given 
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FIG. 18: Optimal EXX OEP potentials for water (cc-pVDZ) obtained from the L-curve analysis 
for each basis set. 

in Table HTT1 An orbital basis if cc-pVDZ is used. We use the a bond length of 1.098A, and 
plot potentials along the bond axis. 





s-exponents 


p-exponents 


N2256 (13s3p) 


2", -4 < n < 8 


2™, -2 < n < 


N2512 (14s5p) 


2™, -4 < n < 9 


2™, -3 < n < 1 


N28192 (18s8p) 


2 n ,-4 < n < 13 


2™, -2 < n < 5 



TABLE I: Exponents used in constructing the potential basis sets used for our N2 calculations. 




cc-pVDZ/H2O02(1 0" 4 ) 

cc-pVDZ/H2O01(10" 4 ) 

cc-pVDZ/cc-pVDZ(0) 



A. LDA 

The results for nitrogen are very much similar to those for water previously discussed. 
All three ET basis sets are capable of producing very oscillatory potentials. The L-curve 
analysis in Fig. (|20l) provides appropriate potentials (Fig. ( !2T|) ). 
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FIG. 19: Spectrum of the approximate Hessian for each basis set. 
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FIG. 20: A regularization 

B. EXX 

The L-curve in Fig. fl23l) and potentials have been presented in the text for this case. 
Included here is the Hessian spectrum, and examples of A-regularized potentials for non- 
optimal A. 
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FIG. 21: Optimal LDA OEP potentials for N2 (cc-pVDZ) obtained from the L-curve analysis for 
each basis set. 
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FIG. 22: Spectrum of the Approximate Hessian for each basis set. 
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FIG. 23: EXX L-curves for N2 in the various potential basis sets. 
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FIG. 24: Optimal EXX OEP potentials for N2 (cc-pVDZ) obtained from the L-curve analysis for 
each basis set. 
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FIG. 25: A-regularized potentials from various points on the L-curve for the calculation cc- 
pVDZ/N2256 

IV. ARGON 
A. EXX 

Here we look at the effects of using a very high quality basis set, the Partridge Uncon- 
tracted 3 basis. For the potential basis set we use the Partridge Uncontracted 3 set itself, 
as well as a 20s ET basis Ar32768 with exponents 2 n , —4 < n < 15. 

The L-curve of Fig. (1271) is interesting as it shows an example where the use of coinciding 
basis sets for orbitals and potential does give rise to nonphysical potentials. This would be 
difficult to predict on the basis of the Hessian spectrum of Fig. [26] alone. 

Also notable from the L-curve is the very well defined EXX OEP energy AE = 5.2 x 10~ 3 
a.u. 



B. Small Potential Basis Sets 

Finally, we consider the situation where we use a basis set that is so small and inflexible 
that it is not able to provide an adequate description of the OEP. The final three figures 
repeat those shown in the paper for argon, however now with the inclusion of the small 
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FIG. 26: Spectrum of the Approximate Hessian for each basis set. 
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FIG. 27: EXX L-curves for Ar (Partridge Uncontracted 3) 



potential basis set Ar64. While this basis set gives significantly different L-curves (Figs. 
(1291) and (ED for LDA and EXX respectively) and potentials (Figs. (ED]) and (E2J)) from 
those larger basis sets, we see that the L-curve analysis still is able to generate a sensible 
potential differing only near the nuclei. 
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FIG. 28: Optimal regularized potentials for each basis set, and an unregularised potential 
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FIG. 29: LDA L-curves for Ar (cc-pVDZ) 



> 
> 



1.6- + 

A 

1.4- * 
1.2- + 



0.8- 



0.6- 



A A + &>0 



24 







-2 
-4 
■S -6 

Q_ 

LU 

D 

_i 

"g -10 

_N 

CO 



CD 



-12 
-14 
-16 
-18 



10" 



cc-pVDZ/Ar64(10~ 7 ) 
cc-pVDZ/Ar1 024(1 0" 14 ) 

■cc-pVDZ/Ar8192(10" 14 ) 
cc-pVDZ/cc-pVDZ(0) 

■ LDA 



10"' 10" 
Radial Coordinate (a.u.) 



10' 



FIG. 30: Optimal regularized LDA potentials for each basis set 
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FIG. 31: EXX L-curves for Ar (cc-pVDZ) 
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FIG. 32: Optimal regularized EXX potentials for each basis set 
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